# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""
import os
cur_path = os.getcwd()
os.chdir( cur_path )

import pandas as pd
import numpy as np
import math
import time
t_start = time.time()

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from scipy.interpolate import griddata


cfx_Dataframe = pd.read_csv(cur_path +'\\data\\cfx_export_0518.csv')    
print('CFX data shape:',cfx_Dataframe.shape)
colnames = cfx_Dataframe.columns.values.tolist()
# print(colnames) 



cfx_Data = cfx_Dataframe.to_numpy()
Data_row=np.size(cfx_Data,0)  #计算 X 的行数
Data_col=np.size(cfx_Data,1)  #计算 X 的列数

cfxpoints = np.zeros((Data_row,3))
for i in range(0,3):
    cfxpoints[:,i] = cfx_Data[:,i]  
    
'check the points'
# highpoint = np.array([]).reshape(0,3)
# medpoint = np.array([]).reshape(0,3)
# lowpoint = np.array([]).reshape(0,3)
# j=0
# for i in range(Data_row):
#     if(cfxpoints[i,2]>0.2 and cfxpoints[i,2]<0.5):
#         j=j+1
#         medpoint = np.concatenate([medpoint,cfxpoints[i,:].reshape(1,3)],axis = 0) 
#     if(cfxpoints[i,2]>0.3):
#         j=j+1
#         highpoint = np.concatenate([highpoint,cfxpoints[i,:].reshape(1,3)],axis = 0) 
#     if(cfxpoints[i,2]<-0.1):
#         j=j+1
#         lowpoint = np.concatenate([lowpoint,cfxpoints[i,:].reshape(1,3)],axis = 0) 

# fig = plt.figure() # 创建一个画布figure，然后在这个画布上加各种元素。
# ax = Axes3D(fig)
# # ax.scatter(cfxpoints[:,0],cfxpoints[:,1],cfxpoints[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# # ax.scatter(medpoint[:,0],medpoint[:,1],medpoint[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# ax.scatter(lowpoint[:,0],lowpoint[:,1],lowpoint[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# ax.scatter(highpoint[:,0],highpoint[:,1],highpoint[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# # ax.scatter(grid_points_valid[:,0],grid_points_valid[:,1],grid_points_valid[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# ax.scatter(grid_points[:,0],grid_points[:,1],grid_points[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# ax.set_xlabel('X') # 画出坐标轴
# ax.set_ylabel('Y')
# ax.set_zlabel('Z')
# ax.set_xlim(-0.35,0.35 )
# ax.set_ylim(-0.35,0.35 )
# ax.set_zlim(-0.2,0.5)
# plt.show()


'set grid'

x_size = (110) *10**-3 #尺寸，单位m
y_size = (110) *10**-3
z_size = (50)  *10**-3

x_coord_start = -0.5 * x_size
y_coord_start = -0.5 * y_size
z_coord_start = 0

grid_dist = 1 *10**-3 #每个正方体cell的边长，单位m


x_num = math.ceil(x_size / grid_dist)
y_num = math.ceil(y_size / grid_dist)
z_num = math.ceil(z_size / grid_dist)
# grid_points = np.zeros((x_num*y_num*z_num,3),dtype = float)


r2 = (0.5 * x_size)**2 #半径平方

'initialize the point in geometry flow'
grid_points = np.zeros((x_num*y_num*z_num,3),dtype = float)
grid_position = np.zeros((x_num*y_num*z_num,3),dtype = int)
grid_num = np.shape(grid_position)[0]
point_index = -1
for m in range(x_num):
    for n in range(y_num):
        for k in range(z_num):
            check = (x_coord_start + m*grid_dist)**2 + (y_coord_start + n*grid_dist)**2
            if(check < r2):
                point_index = point_index + 1
                #存储xyz坐标
                grid_points[point_index] = [m*grid_dist + x_coord_start,    \
                                            n*grid_dist + y_coord_start,    \
                                            k*grid_dist + z_coord_start ]
                #存储mnk值    
                grid_position[point_index]=[m,n,k]
                
#只截取需要的点
grid_points = grid_points[0:point_index]
grid_position = grid_position[0:point_index]


#测试abs(velocity)插值
# count = 0     
# coldcount = 0 
# wallcount = 0
"""
0 1 2   
X Y Z

22-24 25-27 28-30 31-33
u     v     w     abs(Velosity)

3-5 12-14 6-8     9-11  15 
CH4 O2    Density Visco Tem  

16-18
Turbulence Eddy Dissipation

19-21
Turbulence Kinetic Energy

"""




'unpack varibles'
cfx_velo = cfx_Data[:,31]

cfx_V_ABS    = cfx_Data[:,[31,33]]
cfx_U       = cfx_Data[:,[22,24]]
cfx_V       = cfx_Data[:,[25,27]]
cfx_W       = cfx_Data[:,[28,30]]
cfx_CH4     = cfx_Data[:,[3,5]]
cfx_O2      = cfx_Data[:,[12,14]]
cfx_DENS    = cfx_Data[:,6]
cfx_VISC    = cfx_Data[:,9]
cfx_TED     = cfx_Data[:,16]#Turbulence Eddy Dissipation
cfx_TKE     = cfx_Data[:,19]#Turbulence Kinetic Energy

# del(cfx_Data)




'interpolate'
grid_velocity = griddata(cfxpoints, cfx_velo, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
# del(cfx_velo)

grid_V_ABS = griddata(cfxpoints, cfx_V_ABS, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_U = griddata(cfxpoints, cfx_U, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_V = griddata(cfxpoints, cfx_V, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_W = griddata(cfxpoints, cfx_W, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_CH4 = griddata(cfxpoints, cfx_CH4, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_O2 = griddata(cfxpoints, cfx_O2, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_DENS = griddata(cfxpoints, cfx_DENS, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_VISC = griddata(cfxpoints, cfx_VISC, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_TED = griddata(cfxpoints, cfx_TED, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
grid_TKE = griddata(cfxpoints, cfx_TKE, (grid_points[:,0],grid_points[:,1],grid_points[:,2]), method='linear')
   


'output'
delete_=[]
for i in range(np.shape(grid_velocity)[0]):
    if(np.isnan(grid_velocity[i])):
        delete_.append(i)

# grid_velocity_valid = np.delete(grid_velocity,np.array(delete_),axis=0)
# np.save('data\\grid_velocity_valid',grid_velocity_valid)
if(len(delete_) != 0):
    grid_points_valid = np.delete(grid_points,np.array(delete_),axis=0)
    grid_position_valid = np.delete(grid_position,np.array(delete_),axis=0)
    grid_V_ABS_valid = np.delete(grid_V_ABS,np.array(delete_),axis=0)
    grid_U_valid = np.delete(grid_U,np.array(delete_),axis=0)
    grid_V_valid = np.delete(grid_V,np.array(delete_),axis=0)
    grid_W_valid = np.delete(grid_W,np.array(delete_),axis=0)
    grid_CH4_valid = np.delete(grid_CH4,np.array(delete_),axis=0)
    grid_O2_valid = np.delete(grid_O2,np.array(delete_),axis=0)
    grid_DENS_valid = np.delete(grid_DENS,np.array(delete_),axis=0)
    grid_VISC_valid = np.delete(grid_VISC,np.array(delete_),axis=0)
    grid_TED_valid = np.delete(grid_TED,np.array(delete_),axis=0)
    grid_TKE_valid = np.delete(grid_TKE,np.array(delete_),axis=0)
else:
    grid_points_valid = grid_points
    grid_position_valid = grid_position
    grid_V_ABS_valid = grid_V_ABS
    grid_U_valid = grid_U
    grid_V_valid = grid_V
    grid_W_valid = grid_W
    grid_CH4_valid = grid_CH4
    grid_O2_valid = grid_O2
    grid_DENS_valid = grid_DENS
    grid_VISC_valid = grid_VISC
    grid_TED_valid = grid_TED
    grid_TKE_valid = grid_TKE

np.save('data\\grid_points_valid',grid_points_valid)
np.save('data\\grid_position_valid',grid_position_valid)
np.save('data\\grid_V_ABS_valid',grid_V_ABS_valid)
np.save('data\\grid_U_valid',grid_U_valid)
np.save('data\\grid_V_valid',grid_V_valid)
np.save('data\\grid_W_valid',grid_W_valid)
np.save('data\\grid_CH4_valid',grid_CH4_valid)
np.save('data\\grid_O2_valid',grid_O2_valid)
np.save('data\\grid_DENS_valid',grid_DENS_valid)
np.save('data\\grid_VISC_valid',grid_VISC_valid)
np.save('data\\grid_TED_valid',grid_TED_valid)
np.save('data\\grid_TKE_valid',grid_TKE_valid)

grid_num = np.shape(grid_points_valid)[0]
print('valid points: %d' % grid_num)
t_end = time.time()
# print('fitting program ellapsed time: %.2f (s)'%(t_end-t_start))
print('ellapsed time:','%dm'%math.floor((int(t_end)-int(t_start))/60),'%ds'%int((t_end-t_start)%60))
# CFX data shape: (141084, 34)
# fitting program ellapsed time: 71.80 (s)

colname = ['x','y','z']
cell_dataframe = pd.DataFrame(data = grid_points_valid,columns = colname)
cell_dataframe.to_csv('data\\grid_points_valid.csv',index=None)



# fig = plt.figure() # 创建一个画布figure，然后在这个画布上加各种元素。
# ax = Axes3D(fig)
# ax.scatter(cfxpoints[:,0],cfxpoints[:,1],cfxpoints[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# # ax.scatter(medpoint[:,0],medpoint[:,1],medpoint[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# # ax.scatter(lowpoint[:,0],lowpoint[:,1],lowpoint[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# # ax.scatter(highpoint[:,0],highpoint[:,1],highpoint[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# # ax.scatter(grid_points_valid[:,0],grid_points_valid[:,1],grid_points_valid[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# ax.scatter(grid_points[:,0],grid_points[:,1],grid_points[:,2]) # 画出(xs1,ys1,zs1)的散点图。
# ax.set_xlabel('X') # 画出坐标轴
# ax.set_ylabel('Y')
# ax.set_zlabel('Z')
# ax.set_xlim(-0.11,0.11 )
# ax.set_ylim(-0.13,0.13)
# ax.set_zlim(-0.2,0.5)
# plt.show()




# 'test'
# def func(x, y, z):
#     x=x*100
#     y=y*100
#     z=z*100
#     # return (x*(1-x)*np.cos(4*np.pi*x) * (np.sin(4*np.pi*y**2)**2)*z)
#     return 100*x*y+200*z +100*x*y*z
    

# def rmse(predictions, targets):
#     if(np.shape(predictions) != np.shape(targets)):
#         print('rmse input shape error')
#         return np.nan
        
#     for i in range(np.size(predictions,0)):
#         if(np.isnan(predictions[i])):
#             predictions[i] = 0
#             targets[i]=0
#     return np.sqrt(((predictions - targets) ** 2).mean())
#     # return (abs(predictions - targets)).mean()

# cfx_value = func(cfxpoints[:,0], cfxpoints[:,1], cfxpoints[:,2])#实际点的值


# test_calcul_value = griddata(cfxpoints, cfx_value, (grid_points[:,0],grid_points[:,1],grid_points[:,2]),fill_value=np.nan,  method='linear')
# # test_calcul_value = griddata(cfxpoints, cfx_value, grid_points, method='cubic')
# test_actual_value = func(grid_points[:,0],grid_points[:,1],grid_points[:,2])
# rms = rmse(test_calcul_value,test_actual_value)
# std1 = np.std(test_calcul_value) 
# std2 = np.std(test_actual_value)
# mean1 = np.mean(test_calcul_value) 
# mean2 = np.mean(test_actual_value) 
# print('the test of interpolation result \nrmse = ',rms,'\nstdcal=',std1,\
#       '\nstdact=',std2,'\nmeancal=',mean1,'\nmeanact=',mean2)


#反向插值
# test2_calcul_value = griddata(grid_points, test_calcul_value, cfxpoints,fill_value=np.nan,  method='linear')
# test2_actual_value = func(cfxpoints[:,0],cfxpoints[:,1],cfxpoints[:,2])
# rms = rmse(test2_calcul_value,test2_actual_value)
# print('the test of backward interpolation result rms = ',rms)




